home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / src / log.cc < prev    next >
C/C++ Source or Header  |  1997-07-18  |  6KB  |  268 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. #ifdef HAVE_CONFIG_H
  24. #include <config.h>
  25. #endif
  26.  
  27. #include "EIG.h"
  28.  
  29. #include "defun-dld.h"
  30. #include "error.h"
  31. #include "gripes.h"
  32. #include "help.h"
  33. #include "oct-obj.h"
  34. #include "utils.h"
  35.  
  36. // XXX FIXME XXX -- the next two functions should really be just
  37. // one...
  38.  
  39. DEFUN_DLD (logm, args, ,
  40.   "logm (X): matrix logarithm")
  41. {
  42.   octave_value_list retval;
  43.  
  44.   int nargin = args.length ();
  45.  
  46.   if (nargin != 1)
  47.     {
  48.       print_usage ("logm");
  49.       return retval;
  50.     }
  51.  
  52.   octave_value arg = args(0);
  53.  
  54.   int arg_is_empty = empty_arg ("logm", arg.rows (), arg.columns ());
  55.  
  56.   if (arg_is_empty < 0)
  57.     return retval;
  58.   else if (arg_is_empty > 0)
  59.     return Matrix ();
  60.  
  61.   if (arg.is_real_scalar ())
  62.     {
  63.       double d = arg.double_value ();
  64.       if (d > 0.0)
  65.     retval(0) = log (d);
  66.       else
  67.     {
  68.       Complex dtmp (d);
  69.       retval(0) = log (dtmp);
  70.     }
  71.     }
  72.   else if (arg.is_complex_scalar ())
  73.     {
  74.       Complex c = arg.complex_value ();
  75.       retval(0) = log (c);
  76.     }
  77.   else if (arg.is_real_type ())
  78.     {
  79.       Matrix m = arg.matrix_value ();
  80.  
  81.       if (! error_state)
  82.     {
  83.       int nr = m.rows ();
  84.       int nc = m.columns ();
  85.  
  86.       if (nr == 0 || nc == 0 || nr != nc)
  87.         gripe_square_matrix_required ("logm");
  88.       else
  89.         {
  90.           EIG m_eig (m);
  91.           ComplexColumnVector lambda (m_eig.eigenvalues ());
  92.           ComplexMatrix Q (m_eig.eigenvectors ());
  93.  
  94.           for (int i = 0; i < nr; i++)
  95.         {
  96.           Complex elt = lambda (i);
  97.           if (imag (elt) == 0.0 && real (elt) > 0.0)
  98.             lambda (i) = log (real (elt));
  99.           else
  100.             lambda (i) = log (elt);
  101.         }
  102.  
  103.           ComplexDiagMatrix D (lambda);
  104.           ComplexMatrix result = Q * D * Q.inverse ();
  105.  
  106.           retval(0) = result;
  107.         }
  108.     }
  109.     }
  110.   else if (arg.is_complex_type ())
  111.     {
  112.       ComplexMatrix m = arg.complex_matrix_value ();
  113.  
  114.       if (! error_state)
  115.     {
  116.       int nr = m.rows ();
  117.       int nc = m.columns ();
  118.  
  119.       if (nr == 0 || nc == 0 || nr != nc)
  120.         gripe_square_matrix_required ("logm");
  121.       else
  122.         {
  123.           EIG m_eig (m);
  124.           ComplexColumnVector lambda (m_eig.eigenvalues ());
  125.           ComplexMatrix Q (m_eig.eigenvectors ());
  126.  
  127.           for (int i = 0; i < nr; i++)
  128.         {
  129.           Complex elt = lambda (i);
  130.           if (imag (elt) == 0.0 && real (elt) > 0.0)
  131.             lambda (i) = log (real (elt));
  132.           else
  133.             lambda (i) = log (elt);
  134.         }
  135.  
  136.           ComplexDiagMatrix D (lambda);
  137.           ComplexMatrix result = Q * D * Q.inverse ();
  138.  
  139.           retval(0) = result;
  140.         }
  141.     }
  142.     }
  143.   else
  144.     {
  145.       gripe_wrong_type_arg ("logm", arg);
  146.     }
  147.  
  148.   return retval;
  149. }
  150.  
  151. DEFUN_DLD (sqrtm, args, ,
  152.  "sqrtm (X): matrix sqrt")
  153. {
  154.   octave_value_list retval;
  155.  
  156.   int nargin = args.length ();
  157.  
  158.   if (nargin != 1)
  159.     {
  160.       print_usage ("sqrtm");
  161.       return retval;
  162.     }
  163.  
  164.   octave_value arg = args(0);
  165.  
  166.   int arg_is_empty = empty_arg ("sqrtm", arg.rows (), arg.columns ());
  167.  
  168.   if (arg_is_empty < 0)
  169.     return retval;
  170.   else if (arg_is_empty > 0)
  171.     return Matrix ();
  172.  
  173.   if (arg.is_real_scalar ())
  174.     {
  175.       double d = arg.double_value ();
  176.       if (d > 0.0)
  177.     retval(0) = sqrt (d);
  178.       else
  179.     {
  180.       Complex dtmp (d);
  181.       retval(0) = sqrt (dtmp);
  182.     }
  183.     }
  184.   else if (arg.is_complex_scalar ())
  185.     {
  186.       Complex c = arg.complex_value ();
  187.       retval(0) = sqrt (c);
  188.     }
  189.   else if (arg.is_real_type ())
  190.     {
  191.       Matrix m = arg.matrix_value ();
  192.  
  193.       if (! error_state)
  194.     {
  195.       int nr = m.rows ();
  196.       int nc = m.columns ();
  197.  
  198.       if (nr == 0 || nc == 0 || nr != nc)
  199.         gripe_square_matrix_required ("sqrtm");
  200.       else
  201.         {
  202.           EIG m_eig (m);
  203.           ComplexColumnVector lambda (m_eig.eigenvalues ());
  204.           ComplexMatrix Q (m_eig.eigenvectors ());
  205.  
  206.           for (int i = 0; i < nr; i++)
  207.         {
  208.           Complex elt = lambda (i);
  209.           if (imag (elt) == 0.0 && real (elt) > 0.0)
  210.             lambda (i) = sqrt (real (elt));
  211.           else
  212.             lambda (i) = sqrt (elt);
  213.         }
  214.  
  215.           ComplexDiagMatrix D (lambda);
  216.           ComplexMatrix result = Q * D * Q.inverse ();
  217.  
  218.           retval(0) = result;
  219.         }
  220.     }
  221.     }
  222.   else if (arg.is_complex_type ())
  223.     {
  224.       ComplexMatrix m = arg.complex_matrix_value ();
  225.  
  226.       if (! error_state)
  227.     {
  228.       int nr = m.rows ();
  229.       int nc = m.columns ();
  230.  
  231.       if (nr == 0 || nc == 0 || nr != nc)
  232.         gripe_square_matrix_required ("sqrtm");
  233.       else
  234.         {
  235.           EIG m_eig (m);
  236.           ComplexColumnVector lambda (m_eig.eigenvalues ());
  237.           ComplexMatrix Q (m_eig.eigenvectors ());
  238.  
  239.           for (int i = 0; i < nr; i++)
  240.         {
  241.           Complex elt = lambda (i);
  242.           if (imag (elt) == 0.0 && real (elt) > 0.0)
  243.             lambda (i) = sqrt (real (elt));
  244.           else
  245.             lambda (i) = sqrt (elt);
  246.         }
  247.  
  248.           ComplexDiagMatrix D (lambda);
  249.           ComplexMatrix result = Q * D * Q.inverse ();
  250.  
  251.           retval(0) = result;
  252.         }
  253.     }
  254.     }
  255.   else
  256.     {
  257.       gripe_wrong_type_arg ("sqrtm", arg);
  258.     }
  259.  
  260.   return retval;
  261. }
  262.  
  263. /*
  264. ;;; Local Variables: ***
  265. ;;; mode: C++ ***
  266. ;;; End: ***
  267. */
  268.